
clear all;
clc;
load 97_09_final_Q.mat 
% load TranProb_97_09_final_Q.mat
% load TranProb_std_97_09_final_Q.mat
% load InvCDF_97_09_final_Q.mat
freqj=4;% 4 quarters
year_list=[1997:2009]';
year_list_text=num2str(year_list);
quarter_list=[1:4]';
quarter_list_text=num2str(quarter_list);
Rate_list=[0,1:0.5:5];
prob=zeros(13,freqj,10,12);% the last column will be 1, so we ignore the 13th column and set it as 12 columns
prob_std=zeros(13,freqj,10,12); % the 12th column is useless, we define 12 column for coding conveinence.
count=zeros(1,1,10,13);
avgprob=zeros(1,1,10,12); % the last column will be 1, so we ignore the 13th column and set it as 12 columns
stdprob=zeros(1,1,9,10);  % We delete rating 0 in rows and columns and delete withdrawn in columns since we set it as zero bcs of standadization
bin=zeros(9,10); % Rating 0 is deleted so there are 9 lines. Rating 0, Rating 1, withdrawn are deleted so 9 columns left Zeros(yrs, qtrs, ratings, columns)%1-10 column is rating for 0-5, column 11 is default, column 12 is dwithdrawn, column 13 is total
default_prob=zeros(length(year_list)*freqj,3);
InvCDF=zeros(length(year_list)*freqj,3);
Z=zeros(length(year_list)*freqj,3);
PZ=zeros(1,1,9,10);
I=ones(1,1,9,10);
filename='optimization_insam9807_07q4_test1.xls';
% filename2='transition prob.xls';
% filename3='transition prob_std.xls';
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%To calculate quarterly transition prob   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% to calculate probability and standard prob of each cell for first 12 lines.
for i=1:13
for j=1:freqj
    for ii=1:length(Rate_list)
        for jj=1:12
            prob(i,j,ii,jj) = final_result(i,j,ii,jj)/final_result(i,j,ii,13);            
        end
    end 
end 
end

for i=1:13
    for j=1:freqj
        for ii=1:length(Rate_list)
            for jj=1:12
                prob_std(i,j,ii,jj) = prob(i,j,ii,jj)/(1-prob(i,j,ii,12));
            end
        end
%         xx(:,:)=prob_std(i,j,:,:);
%         sheetlist=strcat('',year_list_text(i,:),'.',quarter_list_text(j,:),'');
%         xlswrite(filename3, xx, sheetlist);
%         clear xx sheetlist
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%To calculate the average transition prob %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%avg calculate time
yrwant=11; %2007
freqjwant=3;
%predicted time
yrp=11;
freqjp=4;

for i=2:yrwant-1 %1998~2006
    for j=1:4
        count(1,1,:,:) = count(1,1,:,:)+ final_result(i,j,:,:);
    end
end
for i=yrwant:yrwant %2007
    for j=1:freqjwant %Q3
        count(1,1,:,:) = count(1,1,:,:)+ final_result(i,j,:,:);
    end
end
xx(:,:)=count(1,1,:,:);
xlswrite(filename, xx,'counts');
clear xx

for ii=1:length(Rate_list)
    for jj=1:12 % the last column for total is deleted since all of them are 1
        avgprob(1,1,ii,jj)=count(1,1,ii,jj)/count(1,1,ii,13);
    end
end
xx(:,:)=avgprob(1,1,:,:);
xlswrite(filename, xx,'avgprob');
clear xx

for ii=2:length(Rate_list)% the first row for rating of 0 is deleted
    for jj=2:11% the first column for rating of 0, the 12th column for withdraw are deleted
        stdprob(1,1,ii-1,jj-1) = avgprob(1,1,ii,jj)/(1-avgprob(1,1,ii,12));
    end
end
xx(:,:)=stdprob(1,1,:,:);
xlswrite(filename, xx,'stdprob');
clear xx

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%To calculate the credit score/bins       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xx(:,:)=stdprob(1,1,:,:);
for ii=1:length(Rate_list)-1
    for jj=2:10
        bin(ii,jj)=norminv(sum(xx(ii,jj:10)),0,1);
    end
end
xlswrite(filename, bin,'bins');
clear xx
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%To calculate default probability  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:12 %start from 1998 to 2008
    for j=1:freqj
        default_prob(4*(i-2)+j,1)=year_list(i,1);
        default_prob(4*(i-2)+j,2)=quarter_list(j,1);
        default_prob(4*(i-2)+j,3)=(sum(final_result(i,j,2:10,11))/sum(final_result(i,j,2:10,13)))/(1-(sum(final_result(i,j,2:10,12))/sum(final_result(i,j,2:10,13))));
           %We exclude the debtors whose initial rating is zero           
     end
end
xlswrite(filename, default_prob, 'default_prob');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%To calculate real credit cycle index Z  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=2:12 %start from 1998 to 2008
    for j=1:freqj
        InvCDF(4*(i-2)+j,1)=year_list(i,:);
        InvCDF(4*(i-2)+j,2)=quarter_list(j,:);
        InvCDF(4*(i-2)+j,3)=norminv(default_prob(4*(i-2)+j,3),0,1);
        Z(4*(i-2)+j,1)=year_list(i,:);
        Z(4*(i-2)+j,2)=quarter_list(j,:);
        Z(4*(i-2)+j,3)=-(InvCDF(4*(i-2)+j,3)-mean(InvCDF(1:4*(i-2)+j,3)))/(std(InvCDF(1:4*(i-2)+j,3))); % std is devided by N-1, unbiased           
     end
end
xlswrite(filename, InvCDF, 'InvCDF');
xlswrite(filename, Z, 'Z');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%To calculate forcasted Zhat        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

yhat=-2.105090871026;%(yp fcst inv cdf)
%y(07q4)
%-2.105090871026
% yp
% -2.1474988
% -1.9843944
% -1.9451506
% -1.8807701
% -2.0175352
% -1.9812658
% -2.1171713
% -2.0554904

miu=mean(InvCDF(1:4*(yrwant-2)+freqjwant,3));
sigma=std(InvCDF(1:4*(yrwant-2)+freqjwant,3));
Zhat=-(yhat-miu)/sigma;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%To minimize the distance to get estimate of w  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%x = fzero(@(x) poly(x, b, c), 0)
N=final_result;
P=prob_std(yrp,freqjp,:,:);
X=bin;
wantedZ=Zhat; 
plrfh=@(w)DIF(w,N,P,X,wantedZ);
[foundW,fval,exitflag,output] = fminbnd(plrfh,-0.5,0.5);
fplot(plrfh,[-0.5 0.5],'*')
PZ = calculatePZ(foundW,X,wantedZ);
xx(:,:)=PZ(1,1,:,:);
xlswrite(filename, xx,'PZ');
sheetlist=strcat('',year_list_text(yrp,:),'.',quarter_list_text(freqjp,:),'');
clear xx

%%%%%%%%%%%%%%%%%%%%%%%%%%
%MAE of Kim99 distance   %
%%%%%%%%%%%%%%%%%%%%%%%%%%

naive1=stdprob;
naive2=prob_std(yrwant,freqjwant,2:10,2:11);
n=repmat(N(yrp,freqjp,2:10,13),1,10);
Preal=P(:,:,2:10,2:11);
e_mod=(n.*(Preal-PZ).^2)./(PZ.*(I-PZ));
e_avg=(n.*(Preal-naive1).^2)./(naive1.*(I-naive1));
e_pre=(n.*(Preal-naive2).^2)./(naive2.*(I-naive2));
perf_mod = mae(e_mod);
perf_naive1 = mae(e_avg);
perf_naive2 = mae(e_pre);


disp('fini');

 
